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Figure 2 : The noise level for an imperfect particle distribution for an ellipsoid kernel relative a spherical 
one. In Fig. 2a the corrected error, Eq. 28, is plotted, and in Fig. 2b the uncorrected, Eq. 29. 
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Abstract 

A general expression for the momentum equation in the Smooth Particle 
Hydrodynamics approximation is derived for an arbitrary shape of the kernel. The 
expression is specialized to an ellipsoidal kernel, and compared with its spherical 
counterpart for various degrees of ellipsoidity. For such an ellipsoidal kernel the contribution 
to the adiabatic acceleration from different parts of the kernel is studied, and how the noise 
level in an interaction varies with the shape of the kernel. The two body interaction 
involving kernels with different shapes is also briefly discussed, as well as the conservation 
problems in the sum over the two body interactions. 


1. Introduction 

In Smooth Particle Hydrodynamics, (SPH), the hydrodynamical quantities are 
calculated for a particle and a number of other particles in its neighbourhood, called the 
neighbours. This neighbourhood is defined by a function called the kernel, with an extension 
in space delimited by the smoothing length, h, which confines the neighbours to within a 
certain volume. All lagrangian character of the particles change its set of neighbours each 
iteration as they move around in space. This makes it easy to model three dimensional 
phenomena and avoid mesh tangling problems. The particle formulation also makes it 
possible to describe large density differences in the same model without wasting computer 
resources on empty volumes. The density gradients may be steep in some cases, with a length 
scale of the order of a smoothing length. This may give problems in for example modelling a 
disk, where higher resolution is desired in the vertical rather than the horizontal direction. 

In discussions about SPH, it is sometimes suggested the use of non-spherical kernels 
could offer an improvement of the method. If the kernel is made non-spherical and is allowed 
to attain a shape that is adaptive to its surrounding, according to some method, it would be 
possible to give the method an adaptive resolution for each particle individually. In the case 
of a disk it is desirable with higher resolution in the vertical direction, and therefore it 
could be suitable using an ellipsoidal kernel. Martel et. al. (1994) presented an 
implementation of non-spherical kernels, called Adaptive SPH, (ASPH), which they in 
Shapiro et. al. (1996) used in applications to of cosmological collapse problems. The 
smoothing length is an anisotropic tensor, H, which describes the shape of the kernel. The 
artificial viscosity is often a problem since it, for example, prevents desirable compression in 
the model. H is therefore also used to trace shocks in the model and reduce the use of 
artificial viscosity. Shapiro et. al. argue that their method gives better results than the 
standard formulation of SPH, particularly when the models involve large anisotropic 
density changes and strong shocks. Further ASPH can give structures with few particles, and 
avoid some of the errouneous artificial viscous heating, that is often a problem in collapse 
simulations. 

Fulbright et. al. (1995) present a similar description of a SPH implementation with 
non-sphererical kernels. They show that spherical kernels cause problems when gas is 
compressed or when a star is torn apart by a black hole, and argue that their spheroidal 
kernel gives better results. 
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Neither Martel et. al. (1994) nor Fulbright et. al. (1995) study the continuum limit of 
the method for a non-sphererical kernel, nor do they discuss possible problems when the 
neighbours are not well distributed. Further Fulbright et. al. (1995) presume that the 
neighbours are distributed highly anisotropically in a compression, without discussing its 
implications. How the neighbours are distributed is essential for the accuracy of the method, 
including the case with a spherical kernel, and will be discussed later. Here it is presumed 
that a linear expansion is a good approximation of the density, according to the assumptions 
in the derivations of the SPH approximations from for example Benz (1990). 

The SPH methodology is briefly described below. More comprehensive introductions to 
the method are found in for example Hernquist & Katz (1989), Monaghan (1992) or Steinmetz 
& Muller (1993). From the momentum equation and standard assumptions in the SPH method, 
an expression for the adiabatic acceleation is derived for the continuous case. This and its 
discrete particle version is used to study variations in the ellipsoidicity of the kernel. 


2. The SPH Method 


In SPH the fluid is modelled by particles that carry necessary physical quantities. 
Other quantities are calculated as a sum over the particles in a volume, described by the 
smoothing length, h. To simulate interactions between particles, the particle mass is 
distributed over this volume, called the kernel. Within 2 h there are a number of other 
neighbours, and h is often varied so that the numbers of neighbours is kept constant. Here this 
is less important, because mainly a continuous fluid will be considered. The smoothing length 
should therefore rather be seen as a measure of a typical length scale in the model. At x k a 
quantity can be approximated by an integration over space: 

(f xk ) = jd 3 X 'f x , k w x , k _ xk , (1) 

where w is the so called kernel. The momentum equation: 
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from for example Benz (1990). In Eq. 3 and hereafter the artificial viscosity is neglected, and 
it is also assumed that dv k /dt = [dv k /dt^. 

To calculate the integrations it is suitable to use spherical symmetry. 


°o n In 

§d 3 x = j dr J"d0 j"d0r 2 sin0, (4) 

ooo 

that is an integration over all space. An equation of state, p = (y-i)np, is necessary as well, 
that relates the pressure, internal energy and density. This is to my knowledge the only 
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equation of state used in SPH. This expression for the pressure is substituted into the 
momentum equation, Eq. 2: 
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In principle hte neighbours may have any distribution, but in the SPH approximation it is 
assumed that a linear expansion is a sufficient approximation for any involved quantity. In 
the SPH approximation of the momentum equation it is therefore assumed that in an 

arbitrary point, x k , 
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The calculation is position independent, and therefore the integration in Eq. (3) is performed 
around the origin, and therefore p * ->p 0 and P > —» Pa, and the prime dropped on x k , and 
the prime dropped on the remaining term. Substitution of p = (y-l)up and a continued 
derivation of the adabatic term in Eq. 3 gives 
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The odd terms are dropped during the calculation under the assumption that the derivative 
of the kernel is odd. 

In all calculations the Monaghan and Lattanzio (1985) kernel is used. 
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0 < r < h 
h < r < 2h , 
r > 2 h 
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where in spherical symmetry the smoothing length is dependent on the angles 9 and 0. 
Further, the derivative of the kernel is used: 


dw dw d(r/h ) 

(1 dr r dh ) 

1 dw 

f £f 1 dll ] 

dw 

dx k ~ d(r/h) dx k 

l/t dx k h 2 dx k J 

d(r/h ) 

[h h 2 d^ ) 

d(r/h ) 
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if h(x k ^j = h[-x k y that is if h is even. In Eq. 7 the odd terms can be neglected, and after 
substitution of Eq. 9 into it gives 
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If the kernel is spherical, dhj dx k =0, then Eq. (11) implies that dv k /dt =-(y-l)C k , 
which satisfies the momentum equation, Eq. (5). Hence Eq. (11) satisfies the momentum 
equation if and only if dh/dx k =0, and any additional term in the integrand is errouneous. 

Fulbright et. al. (1995) do not consider the dh/dx k -terms, and Shapiro et. al. (1996) neglect 
them because, as they argue, these terms are small just as they small with spherical kernels. 
In standard SPH dh/dk ^ is should be small, because otherwise there would be conservation 
problems as well as possible violations of the assumptions that a linear expansion of a 
quantity is a sufficient approximations. Their argument should probably be taken to mean 
that the differences in shapes for neighbouring particles are expected to be small. Here, 
however, the smoothing length describes a non-spherical kernel, that is the dh/dx k -terms 
are not negligable for a sufficently deformed kernel, which invalidates their formulation of 
the SPH-approximation. 

It should be emphazised that for the SPH approximations to be valid iv has to be an 
even function, and hence the calculations will be unphysical unless h is even as well. 
Therefore, ellipsoidal kernels only are considered in the rest of this article. 


3. An Ellipsoidal Kernel 


If the kernel is an ellipsoid, its smoothing length can be written h = ^x'x 1 , as a scalar 
product that is, from the equation of an ellipsoid. 



(13) 


where the smoothing length is a function of the spherical coordinates, h =h(h 0 ,9, 0), where 
h 0 in turn is a function of the density and the number of neighbours. The volume of the 
ellipsoid is by definition 
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which gives the relation XxXyXz = 1 / under the assumtion that the volume is constant when 
h 0 is constant. Here the ellipsoid is assumed to be symmetric in the rt/- plane, that is 


X= yv®. 


(15) 


where % = X x = Xu ar| d co = Xz • This gives a new equation for the smoothing length. 
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The smoothing length is substituted into Eq. 11, which gives 
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and in the z-direction: 


5 



X 


dv z _ 3(y-l)C 2 


K 2 K 


dt 


An 


J d 9 J d 0 sin 9 
o o 


cos 2 9 - 


h 0 sin0cos0(ft)-l/ft) 2 ) \ 

[cosin 2 9 + cos 2 9/co 2 ^ ~ sin0 


COS0 


h 0 (ft) sin 2 9 + cos 2 9/co 2 ) 


- 1/2 


3(/-l)C 2 

47T 


7i 2n 


J d 9 J d <j> sin 9 
o o 


cos 2 9 - 

h 0 sin0cos0(ft)- l/ft) 2 ) \ 

(rosin 2 0 + cos 2 0/ft) 2 )^ ~sin0 



COS0 


h o (ft) sin 2 9 + cos 2 9/co 2 ) 

3 (r-i)c z f 


- 1/2 


\ 


d 9 sin 9 cos" 9 


CO-1 CO 


1 -- 

v ft) sin 2 0 + cos 2 9/or y 


3 (r 1 ) C 2 J d 9sin9cos 2 9 


ft ) 3 -1 


1 -- 

v ft) 3 sin 2 0 + cos 2 9 j 


( 21 ) 


In Fig. 1 dv x /dt and dv z /dt are plotted for log 0.5 < log ft) < log 2.0, where y = 5/3 and C x = 
C z = 1. The figure shows that the acceleration does not only strongly deviate from the 
momentum equation for a moderate deformation of the kernel; a numerical solution gives that 
the acceleration even changes sign at ft) = 0.86 and 0 = 1.36 in the x- and z-direction 
respectively. It can be noted that dv x /dt —> +°° if co —> 0, and asymptotically towards -2.67 if 
co—>°°, while dv z /dt —> -2.67 when o —> 0 and dv , /dt -> +°° if o —> °°. Note that for o = 0 
and co = °° the integrands are not piecewice continuous, and therefore not integrable. 

It is possible to calculate the acceleration for a known shape of the kernel, and relate it 
to its spherical counterpart. Hence it is, at least in principle, possible to compensate for the 
error in the momentum equation as a result of the non-spherical kernels. An example of a 
naive scaling function is 
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which gives a corrected acceleration: 
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There is, however, a problem since S vq will go to infinity for a particular value of co. For a 
perfect distribution of the neighbours, the sum over them would go to zero. There are, 
however, always small deviations in the particles' positions, and therefore the errors become 
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large where S pq is large. Therefore more sophisticated methods to compensate for the 
variation should be constructed. 


4. Anisotropic Contribution to the 
Acceleration 


In a model the distribution of the neighbours is not perfect, which gives an error in the 
calculation of the acceleration. Further, the contribution to the acceleration from different 
parts of the kernel is not isotropic. It is therefore of interest to study how this dependence of 
different parts of the kernel varies for different ellipsoidity. 

Considering ellipsoidal kernels, dividing the integrand in Eq. 21 with 


3(y-l)C* 


h 2 n 

J dr J d (j> 


2 . 7i(y-l)h 3 C k 

r 2 sin# = - - ' ’ 


hi 


sin 9 = - ■ 


7tco 3 (y- l)C k sin# 


(<w 3 sin 2 # + cos 2 #) 


3/2 ' 


(25) 


gives a measure of the contribution to the acceleration in the z-direction per unit volume. This 
gives 
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and from Eq. 22 the corresponding corrected expression is: 
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The noise level due to the imperfect particle distribution relative a spherical kernel 
can now be calculated. In the z-direction for the corrected and uncorrected contribution per unit 
volume, the errors are 
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In Fig. 2 E Z/(0/Corr and E Z/Ca are plotted for log0.5 < log® < log2.0. In Fig. 2a E z ^ corr is 
indefinite at <y = 1.36 where dv z /dt =0, so obviously some other method has to be used to 
compensate for the variation in the acceleration if a non-sphererical kernel is used. 
Furthermore it shows that the error grows relatively rapidly when co becomes small. In Fig. 
2b E z 0) grows fast when co is smaller than one, which shows that when the kernel becomes 
flattened, fewer particles contribute to a major part of the acceleration in the z-direction. The 
x- and y- directions give similar relations when the kernel is elongated. This is, however, a 
problem only if the particle distribution gives a more than negligable error in the 
corresponding spherical kernel. This depends on the local particle distribution and the model 
as such. 


5. The Two Body Interaction 

Now we will study the particle discretization, and two nearby particles that interact 
with each other. From Eqs. 3 and 18 the contribution from particle j on particle i is 
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If h j = h j and 5/; /dx k . = dh /<&*■ . in Eq. 22, it is symmetric regarding z and /, that is 

dv k /dt.. = dv k /dt... If on the other hand h ,■ ^ /i ; - or dh/dx k ; =£ dh/dx k ., the momentum is not 

conserved in the two body interaction. (There are thermal energy conservation problems as 
well, but the enthalpy equation is not considered here.) In the standard formulation, (where 

dh/dx k =0) this is usually solved using the arithmetic mean over the smoothing lengths. 
This does, however, introduce a lack of precision in the physical description of the fluid, 
because the smoothing length is a function of the density described by the environment of the 
particle in question. 

The problem can be exemplified by the sum over j in Eq. 30, where all quantities are 
constant except the smoothing length, that varies as 
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(32) 


From the momentum equation, Eq. 2, the acceleration is zero if the pressure is constant. If the 
smoothing length varies over the kernel, this is not generally true in the SPH approximation. 
Hence, the momentum equation may well be violated even if non-spheric kernels are 
introduced. 

Although the smoothing length is the variable that describes the environment in the 
SPH approximation, it is the neighbours and their distribution around the particles that 
define this neighbourhood. It is assumed that a linear expansion of the physical quantities 
the particles carry is a good approximation. Obviously the distribution is not perfect, which 
induces errors in the SPH approximation as discussed by Selhammar (1997). Although the 
averaging of the smoothing length gives an additional error, the neighbour's distribution in 
itself introduces errors in the environment, and therefore it can be questioned if the exact 
geometrical point where the calculation is performed matters. Variation in the smoothing 
length may, however, give large errors, as Hernquist (1993) noted. 

If the shapes of the particles i and j are different, the situation is more complicated 
and the errors possibly larger. It is not sufficient to average over the smoothing lengths. 
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because Eq. 30 includes its derivate, dh/dx k . To conserve the momentum it is possible to 
average over h and dh / dx k , that is to do the transformation 
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As discussed in Sect. 4, the calculated acceleration must be compensated to satisfy the 
momentum equation. As shown there were difficulties with the suggested form, but it should 
be sufficient for an ideal particle distribution considered here. Averaging over the smoothing 
lengths and its derivative gives 
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The acceleration for particle i as a sum over the neighbours from Eqs. 31-34 becomes 

dv k _ m(r-l)ud kpqr y spq dw ^ 

dt i p “ dx 1 ' j 

where it is assumed that all quantities except the smoothing length are constant. Assuming a 
continuous expression for the smoothing length it is possible to rewrite Eq. 35 as an integral. 
From the momentum equation, Eq. 2, the acceleration should be zero since the pressure is 
constant. This is generally not true in an integral version of Eq. 35, neither in the case with 
varying shape of the kernel, nor with spherical kernels but variable smoothing length. The 
expression is not analysed, partly because the suggested scaling is not acceptable, partly 
because a few examples for certain values are not very enlightening. It suffices to note that 
spherical kernels give errors if the smoothing length varies over the kernel, and that a 
similar problem naturally occurs when a varying shape of the kernel is used. The errors do, 
however, become more severe, and at some degree of kernel deformation, it becomes 
physically meaningless to treat the fluid as a sum of two-particle interactions. 

Another concern is the particle distribution, which according to the SPH- 
approximation is assumed to be isotropic. Using a spherical kernel on a cubic centered initial 
distribution, which for example Katz & Gunn (1991) did, there are directions in which 
particles will be found and others where there are no particles. Then a quantity does not 
satisfy the linear approximation, but additional terms are needed. There are to my 
knowledge no calculations on the validity of the linear approximation in a cubic centered grid 
or other particle distributions, and in a model when the particles move from their initial 
positions. However, the linear expansion can be expected to be a sufficient approximation for 
the cubic centered grid. From the figures in both Fulbright et. al. (1995) and Shapiro et. al. 
(1996) the particles are highly anisotropicly distributed, and piled upon each other in the 
direction of compression. Given the arguments above, it is most certainly not possible, neither 
using a spherical nor ellipsoidal kernel, to get an accurate approximation of a quantity in 
such a particle distribution. 


6. Discussion 
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In SPH the acceleration in a point in space is defined as an integration over the kernel 
around the point in question. In the standard SPH approximation this volume is a sphere 
defined by the smoothing length. Here this relation is derived for an arbitrary shape of the 
kernel, as used by for example Martel et. al. (1995) in the discrete particle approximation. 
The relation includes the derivative of the smoothing length regarding the angles in 
spherical symmetry. To satisfy the momentum equation it is shown that the kernel has to be 
a symmetric even function around the point. This implies that the kernel can be described by 
at the most three parameters, such an ellipsoid, and that the tensor H Martel et. al. uses has 
to be diagonal. 

The kernel is then defined to be symmetric in the xy-plane. Thus, it is decribed by its 
semiaxis in the z-direction, and a smoothing length that gives a spherical kernel of the same 
size as the ellipsoidal one. Calculations for some moderately deformed shapes of the kernel 
show that the extra term with the derivative of the smoothing length, Eq. (19), is large for a 
relatively moderate ellipsoidity, as Fig. 1 shows. 

The errouneous extra term can, in principle, be corrected for. It does however affect the 
acceleration to change sign at rather moderate ellipsoidicities. The suggested 
straightforward corrected expressions therefore have singularities at these points. There 
may, however, be other possible solutions to correct for the deviation in the SPH 
approximation from the momentum equation. 

In the integration different parts of the kernel contribute at different amounts to the 
total acceleration. This angular dependent contribution can be calculated for an ellipsoidal 
kernel and compared with its spherical counterpart. The problem with the comparison is 
that the uncorrected calculation is angular dependent in itself, and that the suggested 
correction is not acceptable. Nevertheless, this is sufficient in order to note that the noise in 
an ellipsoidal kernel due to imperfect particle distribution is higher relative to its spherical 
counterpart. To study this problem in more detail requires knowledge about the implications 
due to the particle distribution in itself, a topic not covered here. 

SPH is often considered as a sum of two-particle interactions. This in contrary to an 
interation, or discretization over the neighbourhood of a particle. In a two-body interaction it 
is possible to average parameters such as the pressure or the smoothing length, to conserve 
properties such as energy, momentum or angular momentum. In an example all quantities are 
held constant except the smoothing length, which as other parameters is approximated with 
a linear expansion. Although not a simulation of the two body interaction situation, it 
exemplifies the problem with varying smoothing length. According to the momentum 
equation the calculation should give a zero acceleration, but is a function of the derivative of 
the smoothing length. This is a result equivalent with the earlier stated conclusion, that the 
kernel has to be an even function around the calculation point. The smoothing length is 
related to the third root of the density, and therefore, (provided the density dstribution is 
well behaved in the neighbourhood) the smoothing length should not vary much. In the case 
with varying ellipsoidities, the derivative of the smoothing length between the two 
particles in the interaction can be considerable. 

Calculations of the two artificial viscous acceleration terms by Monaghan and 
Lattanzio (1985) give different dependencies on the smoothing length. This in turn gives 
different angular dependencies for the contributions to the viscous acceleration. Such a 
calculation is straighforward, although, due to the definition of the artificial viscosity, the 
odd terms can not be neglected in a similar way as in the derivation of Eq. 11. This means that 
the kernel is not an even function in the calculation of the artificial viscosity and introduces a 
non-linear term in the acceleration which could give an non-linear velocity distribution. In 
practice, however, this should probably not influence the calculation as a whole. Meglicki et. 
al. (1993) in their Appendix A calculate an estimation of the Reynold's number using the 
linearity properties of the velocity distribution. They do, however, not take into account that 
half of volume of the kernel by definition does not contribute to the artificial viscosity. 
Therefore they cancel the odd terms in the integration and errouneously exclude the linear 
term. Hence they calculate an overestimation of the Reynolds number based upon the 
quadratic velocity terms, which should not be considered in the linearization of the 
environment of a point in the fluid in the first place. 

Without considering the problem of adapting the shape of the kernel to its 
environment, the consequences of a simple implementation of ellipsoidal kernels are 
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discussed. It is shown that the method faces serious problems that have to be solved before 
any attempt to use it in any model. The theoretical aspects of SPH are not fully investigated, 
and some the difficulties discussed here probably originates from a vague description of the 
SPH method as a whole, a lack of test cases and a general preference of producing an efficient 
code rather than an efficient method. 


References 

Benz, W. 1990, The numerical Modelling of Nonlinear Stellar Pulsations, Problems and 
Prospects, Ed. J. R. Buchler, Kluwer Academic Publishers, Dordrecht, Boston, London 
Fulbright, M. S., Benz, W. & Davies, M. B. 1995, ApJ, 440, 254 
Hernquist, L. & Katz, N. 1989, ApJS 70, 419 
Hernquist, L. 1993, ApJ 404, 717 
Katz, N., Gunn, J. E., ApJ 377, 365 

Martel, H., Shapiro, P. R., Villumsen, J. V. & Kang, H. 1995, Mem. Soc. Astron. Ital. 65:4, 
pl061-1071 

Meglicki, Z., Wickramasinghe, D. & Bicknell, G. V., 1993, MNRAS 264: 691-704 

Monaghan, J. J., 1992, A&AR, 30, 543 

Monaghan, J. J. & Lattanzio, J. C. 1985, A&A 149, 135 

Shapiro, P. R:, Martel, H., Villumsen, J. V. & Owen, J. M., 1996, ApJS 103: 269-330 
Selhammar, M., 1997, In preparation 
Steinmetz, M., Muller, E., 1993, A&A 268, 391 


12 



